# -*- coding: utf-8 -*-
#灯光推小区人口
#遍历小区时获取2016灯光、2016人口和2018灯光数据

import utm
import os
import pyproj
from itertools import islice
import numpy as np
from osgeo import gdal
import pandas as pd
from shapely.geometry import Point
from geopandas import GeoDataFrame
from shapely.ops import transform
from shapely.geometry import Polygon
from sklearn.model_selection import train_test_split        #构造训练集和测试集
from sklearn.linear_model import LinearRegression           #引入线性回归

toponode={}
lat,lon = [],[]
dataSet = []

#栅格处理函数，输入矩形经纬度坐标、栅格数据集路径filename
#返回zone_point（经纬度坐标）
def raster_operation(xmin, ymax, xmax, ymin,filename):
    # 读取部分数据集
    # print('******************')
    # print('开始读取区域栅格...')
    # os.chdir(r'E:\MyWork\0_SD_Local\project\03_cq\data')
    in_ds = gdal.Open(filename)
    in_gt = in_ds.GetGeoTransform()  # 获得输入栅格数据的地理变换（geotransform）

    inv_gt = gdal.InvGeoTransform(in_gt)  # 获得地理变换（geotransform)的逆
    if gdal.VersionInfo()[0] == '1':
        if inv_gt[0] == 1:
            inv_gt = inv_gt[1]
        else:
            raise RuntimeError('Inverse georansform failed')
    elif inv_gt is None:
        raise RuntimeError('Inverse geotransform failed')

    offsets_ul = gdal.ApplyGeoTransform(inv_gt, xmin, ymax)  # 获得左上角和右下角偏移坐标（图像坐标）
    offsets_lr = gdal.ApplyGeoTransform(inv_gt, xmax, ymin)
    off_ulx, off_uly = map(int, offsets_ul)
    off_lrx, off_lry = map(int, offsets_lr)
    # print('左上角图像坐标:')
    # print(off_ulx, off_uly)
    # print('右下角图像坐标:')
    # print(off_lrx, off_lry)

    rows = off_lry - (off_uly - 1)  # 获得要提取区域的行数和列数
    colums = off_lrx - (off_ulx - 1)  # 图上坐标中左上角的x,y坐标最小
    # print('行数：' + str(rows), '列数：' + str(colums))

    in_band = in_ds.GetRasterBand(1)
    data = in_band.ReadAsArray(off_ulx, off_uly, colums, rows)  # 研究区域的栅格数据，Numpy二维数组格式
    # print("区域栅格数组：\n"+str(data))

    # 研究区域的栅格数据处理
    # print('开始栅格数据处理...')
    xyCoors = []
    GRID_CODE = []
    for r in range(off_uly, off_lry + 1):  # 获得区域栅格的点坐标以对应的值
        for c in range(off_ulx, off_lrx + 1):
            xyCoors.append(gdal.ApplyGeoTransform(in_gt, c, r))
            GRID_CODE.append(data[r - off_uly, c - off_ulx])

    # print('坐标点个数：')
    # print(len(xyCoors))
    # print(xyCoors)
    # print('点的值个数：')
    # print(len(GRID_CODE))
    # print(GRID_CODE)

    s1 = pd.Series(xyCoors)  # 创建DataFram
    s1.name = 'geometry'
    s2 = pd.Series(GRID_CODE)
    s2.name = 'GRID_CODE'
    df = pd.DataFrame([s1, s2], index=['geometry', 'GRID_CODE'])
    df = df.T

    geometry = [Point(xy) for xy in df['geometry']]  # DataFram转化为GeoDataFram
    df = df.drop(['geometry'], axis=1)  # 从行或列中删除指定的标签。
    crs = {'init': 'epsg:4326'}
    zone_point = GeoDataFrame(df, crs=crs, geometry=geometry)
    # print('栅格处理已完成!')
    # print('******************')
    return zone_point

#经纬度转UTM，输入（经度，纬度）
def wgs2proj(y_84,x_84):
    (x_coor, y_coor, temp_1, temp_2) = utm.from_latlon(x_84, y_84)
    if not (zone_number == temp_1):                                                 # 存在跨zone的点
        p1 = pyproj.Proj(proj='utm', zone=temp_1, datum='WGS84')
        p2 = pyproj.Proj(proj="utm", zone=zone_number, datum='WGS84')
        x_coor, y_coor = pyproj.transform(p1, p2, x_coor, y_coor)                  #转换带号
    return x_coor, y_coor

#文件夹目录处理
f = open("../files/runinfo_vc6.txt")
next(f)
for d in f:
    d = d
coorDic = str(d)+"/data/coor.txt"
latlon_utmDic = str(d)+"/data/latlon-utm.txt"
zoneDic = str(d)+"/data/zone.txt"
toponodeDic = str(d)+"/data/toponode.txt"
TripPopulationDic = str(d)+"/data/TripPopulation.txt"
mainAreaDic = str(d)+"/cache/mainArea.txt"
autoZoomParaDic = str(d)+"/cache/autoZoomPara.txt"


#用户选择的灯光人口数据文件路径
f = open(autoZoomParaDic)
next(f)
d2 = f.readline()
d3 = f.readline()
d4 = f.readline()
pop_lig_type = d2.split()[1]    #数据类型：灯光or人口
pop_lig_file = d3.split()[1]    #数据路径，包含文件名
city_type = d4.split()[1]       #城市类型：0、1、2、3

#2016年灯光和人口数据路径
china_light_2016_Dic = "../files/china_light_2016.tif"
china_population_2016_Dic = "../files/china_population_2016.tif"


f = open(latlon_utmDic)                   #获得zone和band
d1 = f.readline()
zone_number = int(d1.split()[5])
zone_letter = d1.split()[6]

for line in open(toponodeDic):
    toponode[line.split()[0]]=(float(line.split()[1]),float(line.split()[2]))     #将拓扑点的坐标做成一个字典，键为ID，值为坐标的元组

print("--开始遍历小区--")
sumPop,sumLight = {},{}
whiteLight,blackLight = {},{}
sumLight_new = {}
whiteLight_new,blackLight_new = {},{}
for line in islice(open(zoneDic),1,None):
    temp=line.split()
    p=map(lambda x:toponode[x],temp[1:])                                                   #用MAP函数对所有拓扑点ID找到其对应的UTM坐标
    tz = Polygon(list(p))                                                                  #将连续的坐标点转换成面要素

    q = map(lambda x: toponode[x], temp[1:])
    myarray = np.asarray(list(q))
    coor_y = myarray[:,0]
    coor_x = myarray[:,1]
    codeNum = len(myarray)                                                                  #单个小区边界中点的个数

    # UTM转经纬度
    for i in range(0, codeNum):
        lat_temp, lon_temp = utm.to_latlon(coor_y[i], coor_x[i], zone_number, zone_letter)
        lat.append(lat_temp)
        lon.append(lon_temp)
    xmin, ymax = min(lon), max(lat)                                                         # 矩形左上角经纬度坐标
    xmax, ymin = max(lon), min(lat)                                                         # 矩形右下角经纬度坐标
    # print('小区矩形范围：')
    # print('左上角坐标:')
    # print(xmin, ymax)
    # print('右下角坐标:')
    # print(xmax, ymin)

    # 调用栅格数据处理函数
    zone_point = raster_operation(xmin, ymax, xmax, ymin,china_light_2016_Dic)       #16年灯光
    zone_point_pop = raster_operation(xmin, ymax, xmax, ymin, china_population_2016_Dic)           #16年人口
    zone_point_new = raster_operation(xmin, ymax, xmax, ymin, pop_lig_file)      # 18年灯光

    # 经纬度转UTM
    zone_point.geometry = zone_point.centroid.apply(lambda z: transform(lambda x, y: wgs2proj(x, y), z))
    zone_point_pop.geometry = zone_point_pop.centroid.apply(lambda z: transform(lambda x, y: wgs2proj(x, y), z))
    zone_point_new.geometry = zone_point_new.centroid.apply(lambda z: transform(lambda x, y: wgs2proj(x, y), z))
    # print(tz)
    # print(zone_point)

    inzone=zone_point[zone_point.within(tz)]                                                #判断在小区内部的点
    inzone_pop = zone_point_pop[zone_point_pop.within(tz)]
    inzone_new = zone_point_new[zone_point_new.within(tz)]

    sumPop[temp[0]] = sum(inzone_pop.GRID_CODE)                                             #2016年小区人口数据

    whiteLight[temp[0]] = len(inzone[inzone.GRID_CODE.values > 5])
    blackLight[temp[0]] = len(inzone) - whiteLight[temp[0]]
    whiteLight_new[temp[0]] = len(inzone_new[inzone_new.GRID_CODE.values > 5])
    blackLight_new[temp[0]] = len(inzone_new) - whiteLight_new[temp[0]]

    # print(inzone.empty)
    if (inzone.empty):
        sumLight[temp[0]] = 0
        sumLight_new[temp[0]] = 0
    else:
        sumLight[temp[0]] = sum((inzone.GRID_CODE/max(inzone.GRID_CODE))*100)                   #2016年小区灯光数据,2016和2018数据的量程不同，做归一化处理
        sumLight_new[temp[0]] = sum((inzone_new.GRID_CODE/max(inzone_new.GRID_CODE))*100)       #2018年小区灯光数据

    print(str(temp[0])+"小区")
    # print("2016灯光：")
    # print(sumLight[temp[0]],whiteLight[temp[0]],blackLight[temp[0]])
    # print("2018灯光：")
    # print(sumLight_new[temp[0]],whiteLight_new[temp[0]],blackLight_new[temp[0]])

print("--遍历小区结束--")
print("--开始回归分析--")
s1 = pd.Series(sumLight)
s1.name = 'sumLight'
s2 = pd.Series(whiteLight)
s2.name = 'whiteLight'
s3 = pd.Series(blackLight)
s3.name = 'blackLight'
s4 = pd.Series(sumPop)
s4.name = 'sumPop'
df = pd.DataFrame([s1,s2,s3,s4], index=['sumLight', 'whiteLight', 'blackLight', 'sumPop'])
df = df.T

data = df.loc[:,['sumLight', 'whiteLight', 'blackLight', 'sumPop']]#前三列为灯光数据，第四列为人口数据
feature_cols = ['sumLight', 'whiteLight', 'blackLight']
x = data[feature_cols]
y = data['sumPop']

x_train,x_test,y_train,y_test = train_test_split(x,y,random_state=1)

linreg = LinearRegression()
model = linreg.fit(x_train,y_train)
print("常数项")
d = linreg.intercept_
print(d)                                                    #输出回归方程中常数项
print("回归系数")
abc = linreg.coef_
print(abc)                                                  #输出自变量对应的回归系数

y_pred = linreg.predict(x_test)
sum_mean = 0
for i in range(len(y_pred)):
    sum_mean+=(y_pred[i]-y_test.values[i])**2
sum_error = np.sqrt(sum_mean/len(x_test))                   #测试集行数,为总行数的25%
print("RMSE:",sum_error)                                    #打印均方根误差
print("--回归分析结束--")

# 构建回归模型输入数据：2018年各小区灯光统计数据
s1_new = pd.Series(sumLight_new)
s1_new.name = 'sumLight_new'
s2_new = pd.Series(whiteLight_new)
s2_new.name = 'whiteLight_new'
s3_new = pd.Series(blackLight_new)
s3_new.name = 'blackLight_new'
newLigth = pd.DataFrame([s1_new,s2_new,s3_new], index=['sumLight_new', 'whiteLight_new', 'blackLight_new'])
newLigth = newLigth.T
# print("回归模型输入数据：")
# print(newLigth)

#生成2018年各小区人口统计数据
newPop = linreg.predict(newLigth)
# print("回归模型输出数据：")
# print(newPop)
# print(type(newPop))
a = np.arange(1,len(newPop)+1,1)#输出的小区编号

with open(TripPopulationDic, 'w') as f:
    f.write('小区编号 小区常住人口\n')
    for i in range(len(newPop)):
        f.write(str(a[i])+" "+str(int(newPop[i]))+" "+"\n")
f.close()
print("--TriPopulation.txt已生成--")
